#coding=utf-8
from numpy import*
import matplotlib.pyplot as plt
import time
# import kMeans



#读数据
def loadDataSet(fileName):
    dataMat = []   #创建列表，存储读取的数据
    fr = open(fileName)
    for line in fr.readlines(): #读每一行
        line1=line.strip()     #删头尾空白
        curLine = line1.split('\t') #以\t为分割，返回一个list列表
        fltLine = map(float,curLine)#str 转成  float
        dataMat.append(fltLine)     #将元素添加到列表尾
    return dataMat

#算距离：array或matrix类型的向量
def distEclud(vecA, vecB):                  #两个向量间欧式距离
    return sqrt(sum(power(vecA - vecB, 2))) #la.norm(vecA-vecB)

#初始化聚类中心：随机产生区间范围的值作为新的中心点
def randCent(dataSet, k):
    #特征维度
    n = shape(dataSet)[1]
    #创建聚类中心的矩阵 k x n
    centroids = mat(zeros((k,n)))
    #遍历n维特征
    for j in range(n):
        #第j维特征属性值min   ,1x1矩阵
        minJ = min(dataSet[:,j])
        #区间值max-min，float数值
        rangeJ = float(max(dataSet[:,j]) - minJ)
        #第j维，每次随机生成k个中心
        centroids[:,j] = mat(minJ + rangeJ * random.rand(k,1))
    return centroids

#随机抽取两个样本，初始的中心质点（前提是所有的样本没有重复的）
def randChosenCent(dataSet,k):
    m=shape(dataSet)[0] #样本数
    centroidsIndex=[]
    dataIndex=range(m)#样本的序号
    for i in range(k):
        randIndex=int(random.uniform(0,len(dataIndex)))
        centroidsIndex.append(dataIndex[randIndex])
        #删除已经被抽中的样本
        del(dataIndex[randIndex])
    centroids=dataSet[centroidsIndex]
    return mat(centroids)


#标准的kmeans算法:当所有样本的分配结果不发生改变时，终止迭代
def kMeans(dataSet, k,createCent=randChosenCent, distMeas=distEclud):
    m = shape(dataSet)[0]   #样本总数
    #分配样本到最近的簇：存[簇序号,距离的平方]
    clusterAssment = mat(zeros((m,2)))
    #step1:#初始化聚类中心
    centroids = createCent(dataSet, k)
    clusterChanged = True
    #所有样本分配结果不再改变，迭代终止
    while clusterChanged:
        clusterChanged = False
        #step2:分配到最近的聚类中心对应的簇中
        for i in range(m):
            minDist = inf; minIndex = -1  #对于每个样本，定义最小距离
            for j in range(k):  #计算每个样本与k个中心点距离
                distJI = distMeas(centroids[j,:],dataSet[i,:])
                if distJI < minDist:
                    minDist = distJI; minIndex = j  #获取最小距离,及对应的簇序号
            if clusterAssment[i,0] != minIndex: clusterChanged = True
            clusterAssment[i,:] = minIndex,minDist**2 #分配样本到最近的簇
        # print 'centroids=',centroids
        #step3:更新聚类中心
        for cent in range(k):#样本分配结束后，重新计算聚类中心
            #获取该簇所有的样本点
            ptsInClust = dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]
            #更新聚类中心：mean()对array类型求均值，axis=0沿列方向求均值。
            centroids[cent,:] = mean(ptsInClust, axis=0)
    return centroids, clusterAssment

#标准kmeans:当误差平方和SSE收敛，退出循环
def kMeansSSE(dataSet,k,createCent=randChosenCent, distMeas=distEclud):
    m = shape(dataSet)[0]
    #分配样本到最近的簇：存[簇序号,距离的平方]
    clusterAssment=mat(zeros((m,2)))
    #step1:#初始化聚类中心
    centroids = createCent(dataSet, k)
    sseOld=0
    sseNew=inf
    while(abs(sseNew-sseOld)>0.001):
        sseOld=sseNew
        #step2:将样本分配到最近的质心对应的簇中
        for i in range(m):
            minDist=inf;minIndex=-1
            for j in range(k):
                #计算第i个样本与第j个质心之间的距离
                distJI=distMeas(centroids[j,:],dataSet[i,:])
                #获取到第i样本最近的质心的距离,及对应簇序号
                if distJI<minDist:
                    minDist=distJI;minIndex=j
            clusterAssment[i,:]=minIndex,minDist**2 #分配样本到最近的簇
        # print 'centroids=',centroids
        sseNew=sum(clusterAssment[:,1])
        # print 'sse=',sseNew
        #step3:更新聚类中心
        for cent in range(k):#样本分配结束后，重新计算聚类中心
            ptsInClust=dataSet[nonzero(clusterAssment[:,0].A==cent)[0]]
            #按列取平均,mean()对array类型
            centroids[cent,:] = mean(ptsInClust, axis=0)
    return centroids, clusterAssment


#先将所有点作为一个簇，然后将该簇一分为二。之后
#选择误差平方和SSE最小的簇划分为两个簇。以此
#进行下去，直到簇的数目等于用户给定的数目k为止
def biKmeans(dataSet,k,disMeas=distEclud):
    m=shape(dataSet)[0]
    clusterAssment=mat(zeros((m,2)))
    #初始所有样本作为一个簇，仅一个中心点
    centroid=mean(dataSet,axis=0).tolist()[0]
    centList=[centroid]#存一个list类型的中心点，方便用append添加新的中心点
    #遍历所有样本，计算到中心点的距离
    for j in range(m):
        clusterAssment[j,1]=disMeas(mat(centroid),dataSet[j,:])**2
    #绘图：只有一个簇，簇的中心为所有样本的均值
    #plt.figure(len(centList))
    #datashow(dataSet,len(centList),mat(centList),clusterAssment)
    while(len(centList)<k):
        print'the %d iter:'%len(centList)
        minSSE=inf
        #假设对已有的簇分别进行二分聚类后,计算各SSE,
        #哪个值小，就选哪个簇进行二分聚类，使得簇的个数增加一个
        for i in range(len(centList)):
            #第i个簇中的样本
            CurrClusterData=dataSet[nonzero(clusterAssment[:, 0].A==i)[0],:]
            #若将第i个簇利用标准的kmeans聚成两个簇
            SplitTwoCentroids,splitClustAssment=kMeans(CurrClusterData,2)
            #若将第i个簇划分，被划分成两个子簇的总的SSE
            SplitSSE=sum(splitClustAssment[:,1])
            #若将第i个簇划分，其他没有划分的簇的总的SSE
            NoSplitSSE=sum(clusterAssment[nonzero(clusterAssment[:,0].A!=i)[0],1])
            #若将第i个簇划分，划分后总的SSE
            newSSE=SplitSSE+NoSplitSSE
            # print "SplitSSE, and NoSplitSSE: ",SplitSSE,NoSplitSSE
            print 'the SSE of the cluster %d to be splited: '%i,newSSE
            if newSSE<minSSE:
                bestCentToSplit=i #最终确定要划分的簇
                bestNewCents=SplitTwoCentroids #被划分后得到的两个新中心点
                #相应样本的分配结果0/1，样本到对应簇的距离平方
                bestClustAss=splitClustAssment.copy()#副本，防止一个改动会影响另一个
                minSSE=newSSE
        #*本轮最终的划分结果：样本分配的簇的变化
        #将最佳被划分簇的聚类结果为1的类别，更换类别为len(centList),即上一轮最大簇序号+1
        bestClustAss[nonzero(bestClustAss[:,0].A==1)[0],0]=len(centList)
        #将最佳被划分簇的聚类结果为0的类别，更换类别为bestCentToSplit
        bestClustAss[nonzero(bestClustAss[:,0].A==0)[0],0]=bestCentToSplit
        print 'the bestCentToSplit is: ',bestCentToSplit
        # print 'the len of bestClustAss is: ', len(bestClustAss) #被划分的样本数
        #*本轮最后划分结果：增加一个新的中心点
        centList[bestCentToSplit] = bestNewCents[0,:].tolist()[0]#matrix->list
        centList.append(bestNewCents[1,:].tolist()[0])
        #*本轮聚类结果更新：clusterAssment(类别，SSE)
        clusterAssment[nonzero(clusterAssment[:,0].A == bestCentToSplit)[0],:]=bestClustAss
        #每轮迭代绘图
        #plt.figure(len(centList))
        #datashow(dataSet,len(centList),mat(centList),clusterAssment)
    return mat(centList), clusterAssment


#2维数据聚类效果显示
def datashow(dataSet,k,centroids,clusterAssment):  #二维空间显示聚类结果
    from matplotlib import pyplot as plt
    num,dim=shape(dataSet)  #样本数num ,维数dim

    if dim!=2:
        print 'sorry,the dimension of your dataset is not 2!'
        return 1
    marksamples=['or','ob','og','ok','^r','sb','<g'] #样本图形标记
    if k>len(marksamples):
        print 'sorry,your k is too large,please add length of the marksample!'
        return 1

    #绘所有样本
    for i in range(num):
        markindex=int(clusterAssment[i,0])#矩阵形式转为int值, 簇序号
        #特征维对应坐标轴x,y；样本图形标记及大小
        plt.plot(dataSet[i,0],dataSet[i,1],marksamples[markindex],markersize=6)
    #绘中心点
    markcentroids=['dr','db','dg','dk','^b','sk','<r']#聚类中心图形标记
    for i in range(k):
        plt.plot(centroids[i,0],centroids[i,1],markcentroids[i],markersize=15)

    plt.title('cluster result') #标题
    plt.show()



#*****************kmeans++
if __name__=='__main__':
    datamat=mat(loadDataSet('testSet.txt'))

    k=4

    # #=====================标准kmeans
    # runTime=10
    # for i in range(runTime):
    #     plt.figure(i+1)
    #     #利用标准kmeans聚类(kmeans()函数或者kMeansSSE()函数)
    #     mycentroids,clusterAssment=kMeans(datamat,k)
    #     datashow(datamat,k,mycentroids,clusterAssment)



    # #======================二分kmeans聚类
    runTime=1
    for i in range(runTime):
        plt.figure(i+1)
        #二分kmeans函数内部可以绘图查看每一轮迭代的聚类结果
        #二分kmeans的两种实现方式：biKmeans()和biKmeans1()
        mycentroids, clusterAssment = biKmeans(datamat, k)
        #datashow(datamat,k,mycentroids,clusterAssment)




